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1 Introduction 

In the last decades, there has been substantial interest in simple models for 
electron-phonon (el-ph) interaction in condensed matter. Despite intensive 
theoretical efforts, it was not before the advent of numerical methods in the 
1980's that a thorough understanding on the basis of exact, unbiased results 
was achieved. Although at the present our knowledge of the rather simple cases 
of a single carrier (the polaron problem) or two carriers (the bipolaron problem) 
in Holstein and Frohlich models is fairly complete, this is not true for arbi- 
trary band fillings. There is still a major desire to develop more efficient simu- 
lation techniques to tackle strongly correlated many-polaron models, which 
are expected to describe several aspects of real materials currently under 
investigation, such as quantum dots and quantum wires, high-temperature 
superconductors or colossal-magnetoresistance manganites. 

One of the principle problems in computer simulations of microscopic mod- 
els is the limitation in both system size and parameter values. Whereas the 
former can be overcome for the polaron and the bipolaron problem in some 
cases, it is very difficult to obtain results of similar quality in the many-electron 
case. Moreover, many approaches still suffer from severe restrictions concern- 
ing the parameter regions accessible. For example, interesting materials such 
as the cuprates and manganites are characterized by small but finite phonon 
frequencies — as compared to the electronic hopping integral — and interme- 
diate to strong el-ph interaction. Unfortunately, simulations turn out to be 
most difficult exactly for such parameters, and it is therefore highly desirable 
to improve existing simulation methods. 

In this chapter, we shall mainly review different versions of a recently de- 
veloped quantum Monte Carlo (QMC) method applicable to Holstein-type 
models with one, two or many electrons. The appealing advantages of QMC 
over other numerical methods include the accessibility of rather large sys- 
tems, the exact treatment of bosonic degrees of freedom (i.e., no truncation is 
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necessary), and the possibility to consider finite temperatures to study phase 
transitions. The important new aspect here is the use of canonically trans- 
formed Hamiltonians, which permits the introduction of exact samphng for 
the phonon degrees of freedom, enabUng us to carry out accurate simulations 
in practically all interesting parameter regimes. 

Additionally, based on a generalization of the Lang-Firsov transformation, 
we shall present a simple variational approach to the polaron and the bipolaron 
problem which yields surprisingly accurate results. 

The chapter is organized as follows. In section 2, we present the general 
model Hamiltonian. Section 3 is devoted to a discussion of the Lang-Firsov 
transformation, and section 4 contains the derivation of the variational ap- 
proach. The QMC method is introduced in section 5. Section 6 gives a selection 
of results for the cases of one, two and many electrons. Finally, we summarize 
in section 7. 

2 Model 

In this paper we focus on the extended Holstein-Hubbard model defined by 

H = - t ^ c\^Cja + U ^n^inji -|- V n^nj 

i {ij) 
i i 

Here cj^. creates an electron with spin a at site i, and = njo- with 
nicr = c|(jCj^. The phonon degrees of freedom at site i are described by the mo- 
mentum pi and coordinate (displacement) Xi of a harmonic oscillator. The mi- 
croscopic parameters are the nearest-neighbour (denoted by ()) hopping am- 
plitude t, the on-site (Hubbard-) repulsion [/, the nearest-neighbour Coulomb 
repulsion V, the Einstein phonon frequency wq and the el-ph coupling g'. 

This model neglects both long-range Coulomb and el-ph interaction, which 
is often a suitable approximation for metallic systems due to screening. 
Two simple limiting cases of the Hamiltonian (1) are the Holstein model 
(i7 = y = 0) and the Hubbard model {g' = V = 0). In general, the physics 
of the model (1) is determined by the competition of the various interactions. 
Depending on the choice of parameters and band filling, it describes fascinat- 
ing phenomena such as (bi-)polaron formation, Mott- and Peierls quantum 
phase transitions or superconductivity. As we shall see below, the adiabaticity 
ratio 

a = LUo/t (2) 

permits us to distinguish two physically different regimes, namely the adiabatic 
regime a <1 and the non-adiabatic regime a > 1. 

We further define the dimensionless el-ph coupling parameter A = g'"^ / (luoW) , 
where W = 4tD is the bare bandwidth in D dimensions. Alternatively, A may 
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also be written as A = 2Ep/W, i.e., the ratio of the polaron binding energy in 
the atomic limit t ~ 0, Ep ~ g'^/2(jJo, and half the bare bandwidth. A useful 
constant in the non-adiabatic regime is = Ep/uiQ. We exclusively consider 
hypercubic lattices with linear size N and volume N^, and assume periodic 
boundary conditions in real space. 

3 Lang-Firsov transformation 

The cornerstone of the methods presented here is the canonical (extended) 
Lang-Firsov transformation of the Hamiltonian (1). The original Lang-Firsov 
(LF) transformation [1] has been used extensively to study Holstein-type mod- 
els. A well-known, early approximation is due to Holstein [2], who replaced 
the hopping term by its expectation value in a zcro-phonon state, neglecting 
emission and absorption of phonons during electron transfer. However, this 
approach yields reliable results only in the non-adiabatic strong-coupling (SC) 
limit. For A = oo (or t ~ 0), the LF transformation provides an exact solution 
of the single-site problem [3]. 

Whereas transformed Hamiltonians have been treated numerically before 
[4-6], the first QMC method making use of the LF transformation has been 
proposed in [7]. 

We introduce the extended LF transformation by defining the unitary 
operator 

$=e^ , S = iY^ HjuSj (3) 

with real parameters 7ij, i,j = 1, . . . . ^ as defined in equation (3) has 
the form of a translation operator, and fulfills = Given an electron 
at site i, <1> mediates displacements 7jj of the harmonic oscillators at all sites 
j. Hence, the extended transformation is capable of describing an extended 
phonon cloud, important in the large-polaron or bipolaron regime. We shall 
use this transformation for the variational approach. However, the standard 
(local) LF transformation will be expedient as a basis for unbiased QMC 
simulations, in which the transformed Hamiltonian is treated exactly. 

Operators have to be transformed according to vl = ^A^^ . Defining the 
function /(ry) = e'^^Ae~'^^ we obtain 

f'{v)=e^^[S,A]e-^^, (4) 
where /' = df/dt]. A simple calculation gives 

[S,Cicr] ^ -i^ -in Pi Ci„ , [S, c|„] = i ^ 7,, pi . (5) 
( I 

Substitution in equation (4), integration with respect to r] and setting r] = 1 
results in 
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cL = 4a . = Cia e"' . (6) 

For phonon operators, the relation 

1 = e^Ae-^ = A+[S,A] + [5, A]] + • • • , (7) 

yields 

Xi=Xi + ^ ^ijHj , Pi =Pz- (8) 

3 

Collecting these results, the transformation of the Hamiltonian (1) leads to 
H=-tJ2 cLc,.,e'^'(^"-''^')^^' + ^ ^(^2 ^ ^1) + ^ nMcJolij - g'Sij) 

{ij)<T i ij 

V / V / ^ ' 



He. 



Here the term Hep describes the coupling between electrons and phonons, 
whereas Hee represents an effective el-el interaction. Hamiltonian (9) will be 
the starting point for the variational approach in section 4. 

For QMC simulations, it is more suitable to require that the el-ph terms 
in Hep cancel. This can be achieved by setting 7^ = 'ySij with 



7=,^. (10) 

Wo 



The parameter 7 corresponds to the distortion which minimizes the potential 
energy of the shifted harmonic oscillator Spot = — d'x. This leads us to 
the standard LF transformation 

^o = e^°, So = h^niPi, (11) 



and the familiar results for the transformed operators 

4 = 4e'^?S Ci, = ci„e-'^P* (12) 

and 



Xi=Xi+ irii , Pi=Pi- (13) 

In contrast to the non-local transformation (3), only the oscillator at the 
site of the electron is displaced. The transformed Hamiltonian reads 
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{ij)a i 

^ T' ' ' 

+ {\U-E^)Y,nl + vY,nini-\uY,ni. (14) 

i (ij) i 



As we shall discuss in detail in section 5, the difficulties encountered in 
QMC simulations of the original Hamiltonian (1) are to a certain extent re- 
lated to (bi-)polaron effects, i.e., to the dynamic formation of spatially rather 
localized lattice distortions which surround the charge carriers and follow their 
motion in the lattice. 

For a single electron, the aforementioned Holstein- Lang-Firsov (HLF) ap- 
proximation [2] becomes exact in the non-adiabatic SC or small-polaron limit, 
and agrees qualitatively with exact results also in the intermediate-coupling 
(IC) regime [8]. Although it overestimates the shift 7 of the equilibrium posi- 
tion of the oscillator in the presence of an electron, and does not reproduce the 
retardation effects when the electron hops onto a previously unoccupied site, 
the approximation mediates the crucial impact of el-ph interaction on the lat- 
tice. Consequently, the transformed Hamiltonian (14) can be expected to be a 
good starting point for QMC simulations, which then merely need to account 
for the rather small fluctuations around the (shifted or unshifted) equilibrium 
positions. In principle, it would also be possible to develop a QMC algorithm 
based on the Hamiltonian (9) — the basis of our variational approach — with 
the parameters jij determined variationally, but the local LF transformation 
proves to be sufficient. 

The Hamiltonian (14) does no longer contain a term coupling the electron 
density n and the lattice displacement x. By contrast, the extended transfor- 
mation does not eliminate the interaction term completely. On top of that, 
the hopping term involves all phonon momenta Pi as well as the parameters 
jij, and the el-el interaction becomes long ranged [cf equation (9)]. 

For spin dependent carriers with nf ^rii, the interaction term Hcc contains 
a Hubbard-like attractive interaction. Whereas the latter can be treated ex- 
actly in the case of two electrons (section 5.1), the many-electron case requires 
the introduction of auxiliary fields which complicate simulations. However, no 
such difficulties arise for the spinless Holstein model considered in section 6. 

4 Variational approach 

For simplicity, we shall restrict the following derivation to one dimension; an 
extension to D > 1 is straight forward. Furthermore, we only consider finite 
clusters with periodic boundary conditions, altlioiigli infinite systems may also 
be treated. The results of this section have originally been presented in [7, 9]. 
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4.1 One electron 

As noted before, the simple variational method presented here is based on 
the extended transformation (3), leading to the Hamiltonian (9). We treat 
the 7ij as variational parameters which arc determined by minimizing the 
ground-state energy in a zero-phonon basis in which {Hep) = 0. 

For systems with translation invariance the displacement fields satisfy the 
condition 7,,, = Together with X^^nj = 1 for a single electron we get 

^ee = lEi7f-fl'70. 

The eigenvalue problem of the transformed Hamiltonian (9) is solved by 
making the following ansatz for the one-electron basis states 

|lO=cL|0)®ni4'')), « = l,...,iv|, (15) 

where \4>l^^) denotes the ground state of the harmonic oscillator at site ly. The 
non-zero matrix elements of the hopping term are 

{l\ iJkin \l') = -tS^w) n('/'o'^ |e'(^--^"^)?l4^)) 

= -t,5^„,^e-3 , (16) 

where 4'{x) denotes the real-space wavcfunction of the harmonic-oscillator 
ground state. The Kronecker symbol Si^w) forces I and V to represent nearest- 
neighbor sites. A simple calculation gives for the other terms in equation (9) 

(Z|iJph|0=<5«'Y, (/|i?ep|O=0, {l\Hee\l')=5u'(^Y.^^-9'l)j. 

(17) 

In the zero-phonon subspace spanned by the basis states (15), the eigenstates 
of Hamiltonian (9) with momentum k are 

\'^k) = cl\0)®\{\<t>t'^) (18) 

with eigenvalues 

£;(fc)=SHn + ^ + ^E^?-5'70 (19) 

I 

and the kinetic energy 

£;ki„ = -tYj e"=*e-3 1:.(7.-7.+^)' . (20) 

5=±1 
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Defining the Fourier transform 

% = ;^E«"'^' (21) 

and using (7; e K) 

^ 7,7^+5 = J2 %7-ge'«^ = E ^9 cos 5^ ' (22) 

1/ q 9 

we may write 

Ei^in = -i^e'**e-5 2:,(i-=o«9«)7^ =£Q(fc)e-iE,(i-cos9)7^ (23) 
5 

with the tight-binding dispersion eo(A;) = — 2tcosA;. Hence the ground-state 
energy becomes 

^(^) = -(^) + ^ + ?E^^4e%- (24) 

q ^ q 

The variational parameters 7^ are determined by requiring 

3E 0' I 

^ = -7pe(fc)(l - cosp) a;o7p - ■y= = , (25) 

so that the optimal values 7p can be obtained from 

y/Nuo + e{k){l-cosp) ■ ^^^^ 

Since £(fc) depends implicitly on the 7^, equation (26) has to be solved sclf- 
consistently. It has the typical form of the random-phase approximation since 
a variational ansatz for the untransformed Hamiltonian may be written as 

#t|V;,) = 1 ^e'^^^-^e-'^'-^^'^^' |0)®ni<^o''^)' (27) 
i V 

with as defined in equation (3). 

We shall also calculate the quasiparticle spectral weight for momentum 
k = Q, defined as 

Vi^^= (0|?fe=o,^|V-'o) . (28) 



Here \'tpQ) denotes the ground state with one electron of momentum p = Q and 
the oscillators in the ground state |(?!>o). Fourier transformation and the same 
manipulations as in equation (16) lead to 
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i 

= e"3^<i'^?. (29) 

Just as the HLF approximation, the present variational method becomes 
exact in the non-interacting limit (A = 0) and in the non-adiabatic SC limit. 
Furthermore, it yields the correct results both for a = (classical phonons) 
and a = oc', and also gives acciiratc results for large a and finite A, since 
the displacements of the oscillators — only local and generally overestimated 
in the HLF approximation — are determined variationally. 



4.2 Two electrons 

As in the one-electron case, the use of a zero-phonon basis leads to (Hep) = 
Ojind, neglecting the ground-state energy of the oscillators, we also have 
(^ph) = 0. Hence, H = i?kin + -ffee with the transformed hopping term 

^kin = -ieff = ^(^) ^ka^^ka (30) 

{ij}(T ka- 

and e{k) = —2 toff cos(fc). Here the effective hopping 

*eff=^ E e"'^'^"""'"^'^'' (31) 

where rotational invariance has been exploited. For two electrons of opposite 
spin (i.e., nia-rija = for i ^ j) and F = 0, Hee in equation (9) reduces to 

= 2uo - [/ 2 ^ VijUifTlji , Vij = ^Y ^''J^'''' ~ ^'^iJ + l^ijU ■ (32) 
ij I 

The eigenstates of the two-electron problem have the form 

\^k}=Y^A-pA^\0} ^ (33) 
p 

suppressing the phonon component [cf equation (18)], and may be written as 

l^'^) = E e'''^ E 4i4+i, |o) , (34) 

with the Fourier transform 

d = Fd, {F)ip = e'^'f /ViV . (35) 
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The normalization of equation (33) reads 

{^k\^k) =Y.\M^ ■ (36) 
p 

Using equation (33), wc find for the expectation value of -ffkin 

pp' g 

\v ^ ' V ^ ' 

= Y.\~d,\^[e{p)+e{k-p)] 
p 

= -4teff<i^Tfcd. (37) 

In the last step we have introduced = diag[cos(p) + cos(fc — p)] and 
made use of equation (35). 

The expectation value of the interaction term, best computed in the real- 
space representation (34), takes the form 



ij j'j" W 

„t J 



I S"'+iTS"'i"'iT"'jiS""iS"+i'T 1^) 



(2z;o-C/)X|di|2 + l^t;,+,,,|d, 
jl 



= {2vo - U)d}d + 2SVd , (38) 

with the diagonal matrix Vij = SijVi. 

The minimization of the total energy with respect to d yields the eigenvalue 
problem 

(-4ieff Tk + 2V)d= (Eo -2vo + U)d. (39) 

The vector of coefficients d and thereby the ground state are found by min- 
imizing the ground-state energy Eq through variation of the displacement 
fields 7y. Similar to the one-electron case, this procedure takes into account 
displacements of the oscillators not only at the same but also at surround- 
ing sites of the two electrons, and is therefore capable of describing extended 
bipolaron states (see section 6.2). Note that the two-electron problem is diag- 
onalized exactly without phonons (i.e., for A = 0). 



10 



Martin Hohenadler and Wolfgang von der Linden 



5 Quantum Monte Carlo 

In this section, we present an overview of our recently developed QMC algo- 
rithms for Holstein-type models [7,9-11]. 

As mentioned before, in contrast to the variational approach, the QMC 
approaches discussed here, based on the local LF transformation (14) which 
does not contain any free parameters, are unbiased. They yield exact results 
with only statistical errors that can in principle be made arbitrarily small. 

The motivation for the development of improved QMC schemes for Hol- 
stein models stems from the fact that calculations with existing methods often 
suffer from strong autocorrelations, i.e., non-negligible statistical correlations 
between successive MC configurations [7,12]. In fact, autocorrelations may 
render accurate simulations impossible within reasonable computing time. As 
discussed in [7] , the problem becomes particularly noticeable for small phonon 
frequencies and low temperatures. 

Whereas autocorrelations can be avoided to a large extent for one or 
two electrons by integrating out the phonons analytically, no efficient gen- 
eral schemes exist for finite charge-carrier densities (see discussion in [7]). 

In the sequel, we present a general (i.e., applicable for all densities) solu- 
tion for this problem in several steps. First, the effects due to el-ph interaction 
are separated from the free lattice dynamics by means of the LF transforma- 
tion (14). Since the latter contains the crucial impact of the electronic de- 
grees of freedom on the lattice, simulations may be based only on the purely 
phononic part of the resulting action. The fermionic degrees of freedom can 
then be taken into account exactly by reweighting of the probability distri- 
bution. Consequently, we may completely ignore the electronic weights in the 
updating process, and thereby dramatically reduce the computational effort. 
The principal component representation of the phonon coordinates allows ex- 
act sampling and avoids any autocorrelations. 

5.1 Partition function 

We begin by deriving the partition function for the case of a single electron. 
Then we discuss the differences occurring in the cases of two or more carriers. 

One electron 

The partition function is defined as 

Z = Tre-'^" (40) 

with H given by equation (14) and the inverse temperature P = (fceT)"^. 
For a single electron, i?ee = —Ep becomes a constant which needs only to be 
considered in calculating the total energy. 

Using the Suzuki- Trotter decomposition [12], we obtain 
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e-^" « (e-^-^^«ne-^"^?''e-^^^p>')i = , (41) 

where At = (}/L <C 1. Splitting up the trace into a bosonic and a fermionic 
part and inserting L—1 complete sets of oscillator momentum eigenstates we 
find the approximation 

ZL=TrfJ dpidp2 ■ ■ ■ dpL {pi\U \p2) ■ ■ ■ {pl\U \pi) (42) 

with dpr = rii'^K.T"- Es^ch matrix clement can be evaluated by inserting a 
complete set of phonon coordinate eigenstates / da;|a;)(a;|, since all x-integrals 
are of Gaussian form and can easily be carried out. The result is 

(p,|e-^^<|p.+x) = C^°e--^^.(^--^--)\ C = ^f^. (43) 

The normalization factor in front of the exponential has to be taken into 

account in the calculation of the total energy, but cancels when wc measure 
other observables. With the abbreviation Vp = dpidp2 • ■ - dpL the partition 
function finally becomes 

Zl = C'^°'' j Vpw^Wi, (44) 

where 

L 

Wb = e-^^^^ Wi = TviQ, ^ = n e"^^-^^'" • (45) 

r=l 

Here H^^^ corresponds to i?kin with the phonon operators Pi, pj replaced by 
the momenta pi^r, Pj,T on the rth Trotter slice, and its exponential may be 
written as 

e-'^^^iVn = DrKDl , Kjj> = (e^^*''"^) , (£)^),y = %-e'Tf^- , (46) 

^ ^ jj' 

where h^^ is the A*"^ x A^^ tight-binding hopping matrix. To save some com- 
puter time, we employ the checkerboard breakup [13] 

iij) 

Using equation (47), the numerical effort scales as A''^^ instead of N^^ (see 
also section 5.6), but the error due to this additional approximation is of the 
same order Z\t^ as the Trotter error in equation (41). 

According to equation (46), we have the same matrix k for every time slice, 
which is transformed by the diagonal unitary matrices Dj- . The matrix i? can 
be calculated in an efficient way by noting that the transformation matrices 
Dl and Dr+i at time slice r may be combined to a diagonal matrix 
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{D,,r+i\j = 5ye'^(f--+i-f*-) . (48) 

Due to the cyclic invariance of the fermionic trace, Di can be shifted to the 
end of the product, where it combines with d\ to -Di,,i. Hence we can write 

L 

Q=\{KDr,r+l, (49) 

T = l 

with periodic boundary conditions in imaginary time. In the one-electron case, 
the fermionic weight Wf = (n| i? \n) is given by the sum over the diagonal 
elements of the matrix representation of Q in the basis of one-electron states 
(dropping unnecessary spin indices) 

\n) = ci |0) . (50) 

The bosonic action in equation (45) contains only classical variables: 

Sb = Y T.Plr + E (f'*.- - P^'^+lf , (51) 

i,T i,r 

where the indices i = 1, . . . , N'^ and r = 1, . . . , L run over all lattice sites and 
time slices, respectively, and Pi,L+i = Pi,i- It may also be written as 

S^ = J2pJM (52) 

i 

with Pj = (pi,i, . . . ,Pi,L) and a periodic, tridiagonal L x L matrix A with 
non-zero elements 

Since Zl is a trace, it follows that {A)i^l = (^)l,i = —{2u)qAt'^)~^. 
Two electrons 

In contrast to [9] , here we also take into account nearest-neighbour Coulomb 
repulsion V. For two electrons, the Hamiltonian (14) simplifies to 

H = #kin+^ph+-ffee-2i;p , #ee = (?7-2^p) ^ ni|ni4 +y ^ " 

i (ij) 

Again, the constant shift can be neglected in the QMC simulation, but in 
contrast to the single-electron case, we have a non-trivial interaction term. 
The Suzuki- Trotter decomposition yields 



(55) 
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Using the same steps as above we obtain 

L ^ ^ 

w,b = e-^^^\ Wf = Trfn, 12 = [] c-^^^^r^e"^^^- , (56) 

with 5b given by equation (51). 

As pointed out in [9] , the numerical effort for two electrons increases sub- 
stantially in higher dimensions. Therefore, we restrict ourselves to D = 1. 

Previously, we only considered the case of two electrons of opposite spin 
(forming a singlet) [9] . Here we shall also present results for the triplet state. 

Singlet 

In the singlet case we choose the two-electron basis states 

^^\l) = \i,j) = cl^c]^\0) , i,j = l,...,N} , (57) 

where we have used a combined index 1 = 1,..., N'^. The tight-binding hop- 
ping matrix, denoted as k, has dimension x 7V^, and the corresponding 
exponential in equation (56) can again be written as e-^^-^^kL = D^wDt [cf 
equation (45)], where 

{Dr)u' = ^„,e^T(f--+f^-) (58) 

is diagonal in the basis (57). 

The remaining contribution to f2 comes from the effective el-el interaction 
term Hee in terms of the sparse matrix 

(V)«' = ^{Sik e-^-(^-^^-)^^0;fc(e-^^^'<«>)fe;' . (59) 

k 

The momenta p merely enter the diagonal matrix D: the A'^^ x matrices 
V and K are fixed throughout the entire MC simulation. Finally, we have 

n = '[[DrKDlV, (60) 

r 

and the fermionic trace can be calculated as the sum over the diagonal ele- 
ments of the matrix ^2 in the basis (57), i.e., 

Trff? = ^(j,i|r2|i,i) . (61) 

ij 

Triplet 

For two electrons with parallel spin we use the basis states 

{\l)=\i,j)=4cl\0) , t = l,...,N,j = i + h...,N}, (62) 

i.e., double occupation of a site is not possible. Since we can further not 
distinguish between the states and the dimension of the electronic 
Hilbert space is reduced from A^^ (singlet case) to N{N — l)/2. Consequently, 
for the same system size, simulations for the triplet case will be much faster. 
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Many-electron case 

The one-electron QMC algorithm can easily be extended to the spinless Hol- 
stein model with many electrons. For the latter, assuming V = 0, the in- 
teraction term in equation (14) reduces to Hee = —Ep^^ni. Therefore, the 
grand-canonical Hamiltonian becomes 

H = H - nY^rii = -tY, clc^e'-'^P' "^^ ) +iJph - {Ep + m) ^ , (63) 

i {ij) i 

^ ' ^ ' 

where /x denotes the chemical potential. For half filling n = 0.5 [N/2 spinless 
fermions on N sites, cf equation (84)], the latter is given by /x = —E-p, whereas 
for n 7^ 0.5, it has to be adjusted to yield the carrier density of interest. 

The approximation to the partition function may again be cast into the 
form of equation (44), with as defined by equations (45) and (51), respec- 
tively. The fermionic weight is given by 

Wi = Trf(BiB2 ■■■Bl), % = e-^^^^ine"^^-^" . (64) 

Following Blankenbecler et al. [14], the fermion degrees of freedom can be 
integrated out exactly leading to 

Wi = det(l + Br-- - Bl) = det(l + Q) , (65) 

where the A^^ x matrix Br is given by 

Br=DrnDlV. (66) 

Here /t and Dt are identical to equation (46), and 

(V)y = 5ij e^-(^P+'') . (67) 

There is a close relation to the one-electron Green function 

Gij = {c.c]) + (4c^) . (68) 

In real space and imaginary time, we have [14, 15] 

Gtj = {c,cl) = {l + f2)7.^ , Gtj = 5,j - G«. = (12 G''),, . (69) 

At this stage, with the above results for the partition function, a QMC 
simulation of the transformed Holstein model would proceed as follows. In 
each MC step, a pair of indices {io,To) on the A''^ x L lattice of phonon 
momenta pi^r is chosen at random. At this site, a change Pi^.To ^ Pia.ra + '^P 
of the phonon configuration is proposed. To decide upon the acceptance of 
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the new configuration using the Metropolis algorithm [12], the corresponding 
weights wi-,W[ and w'y^w'^ have to be calculated. Due to the local updating 
process, the computation of the change of the bosonic weight Aw^, = w'^^/w^ 
is very fast, which is not the case for the fermionic weight Awf = w'f/wf. 
By varying tq sequentially from 1 to L instead of picking random values, the 
calculation of the ratio of the fermionic weights can be reduced to only two 
matrix multiplications. 

It turns out that a local updating as described above does not permit 
efficient simulations for small phonon frequencies or low temperatures. There- 
fore, we shall introduce an alternative global updating in terms of principal 
components in section 5.4. 

5.2 Observables 

Using the transformed Hamiltonian (14), the expectation value 



As a result of the analytic integration over the phonon coordinates x, 
interesting observables such as the correlation function (riiXj) are difficult to 
measure accurately. Other quantities such as the quasiparticle weight, and the 
closely related effective mass [16], can be determined from the one-electron 
Green function at long imaginary times [17], but results for one electron or 
two electrons would not be as accurate as in existing work (e.g., [18-20]). 

The situation is strikingly different in the many-electron case, for which 
many methods fail to produce results of high accuracy for large systems and 
physically relevant parameters. Moreover, other important observables, such 
as the one-electron Green function, can be calculated with our approach. 

One electron 

The electronic kinetic energy is defined as 



(O) =Z-^Trd e-^" = ^"^ Tr O e"'^^ 



(70) 



of an observables O is computed according to 




(71) 



(72) 



{ij) 



Repeating the steps used to derive the partition function, and noting that the 
additional phase factors in equation (72) again lead to the same matrix i? as 
in equation (49), we find 
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= -tZl^ I 'Dpw^,{j\Q\i) (73) 

with the one-electron states (50). Introducing the matrix elements = 
{i\ Q \j) and the expectation value with respect to 

J Vpwx, 

we obtain 

_ Efa) (%)b 

i-kin = -i ^ ,^ , — • (75) 

Here we have anticipated the reweighting discussed in section 5.3. 
The total energy can be obtained from E = —d{\n.Z)/dj3 as 

i 

To compare with other work we subtract the ground-state energy of the 
phonons, -Bo.ph = 

Two electrons 

For two electrons, exploiting spin symmetry, we have 

i^kin = E = -2t E(4tSt«'"^''""^'^) • (77) 



(78) 



A simple calculation gives 

Writing out explicitly the fermionic trace we obtain 

vy 

= E(j,/|/2|i,/) , (79) 

i' 

and the kinetic energy finally becomes 
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iJkin = -2tZl^ j ^^e'^(f-i-P-i) {j,3'\ n . (80) 

^ m r 

In addition to .Bkin, we shall also consider the correlation function 

p{5) =Y,{ni^ni+si) , <5 = 0,l,...,7V/2-l (81) 

i 

depending on the distance 5. We find 

p{5)=Z£^ f Vpwh^{i,i + 5\n\i,i + S) . (82) 

i 

Many-electron case 

The calculation of observables within the formalism presented here is similar 
to the standard determinant QMC method [13-15]. For an equal-time (i.e., 
static) observable O we have 



/rt\ I T^pwbWfTTf{OBi---BL) , . 

i'-'/h = • (83) 



The carrier density 



^ E("^) (84) 



may be calculated from G'' [equation (69)] using (n,) = (G^^). 

Similarly, the modulus of the kinetic energy per site is given by 

^^^^ = ^E((^%)- (85) 

Equal-time two-particle correlation functions such as 

Pi^) = X](^'"''+*) (8^) 

i 

may be calculated in the same way as in [14, 15]. For a given phonon config- 
uration. Wick's Theorem [3] yields 

{ninj)p = {clcic]cj)p 

and {niTij) is then determined by averaging over all phonon configurations. 
The time-dependent one-particle Green function 
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G\k,T) ^ {cI{t)c, 



ki — \^ ^fe 



is related to the momentum- and energy-dependent spectral function 

A{k,uj- fj) = --lm.G^{k,(jO- jj) (89) 



through 



The inversion of the above relation is ill-conditioned and requires the use of 
the maximum entropy method [11, 12,21]. Fourier transformation leads to 



ik,r) = ^Y.^"-^^'~^'^Gtj{r). (91) 



The allowed imaginary times are r; = IAt, with non-negative integers < 
I < L. Within the QMC approach, we have [14, 15] 

G^,.(rO = (G"i3i...i?0,.. (92) 
The one-electron density of states is given by 

N{w - 1^) = --Im G{uj - n) , (93) 

TT 

where G{iv — /z) = (N^)^^ J2k G{k,u) — /z). It may be obtained numerically 
via 

N{T)=GUr), (94) 
and subsequent analytical continuation. 



Suzuki- Trotter error 



The error associated with the approximation made in, e.g., equation (41) 
can be systematically reduced by using smaller values of At. In practice, 
there are two strategies to handle this so-called Suzuki- Trotter error. Owing 
to the usually large numerical effort for QMC simulations, At is often simply 
chosen such that the systematic error is smaller than the statistical errors for 
observablcs. A second, more satisfactory, but also more costly method is to 
run simulations at different values of At, and to exploit the At"^ dependence 
of the results to extrapolate to At = 0. 

For the results in section 6, wc have used a scaling toward At = based 
on typical values At = 0.1, 0.075 and 0.05 to obtain the results for one and 
two electrons. In contrast, for the numerically more demanding calculations 
of dynamic properties in the many-electron case, At = 0.1 has bc;en chosen. 
This is justified by the uncertainties in the analytical continuation. 
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5.3 Reweighting 

As pointed out at the end of section 5.1, the calculation of the change of the 
fermionic weight W{ represents the most time-consuming part of the updating 
process. Consequently, it would be highly desirable to avoid the evaluation of 
Wf. This may be achieved by using only the bosonic weight w^, in the updating, 
and treating wt as part of the observables. For the expectation value of an 
observable O, such a reweighting requires calculation of 

(0) = ^, (95) 

where the subscript "b" indicates that the average is computed based on Wh 
only [cf equation (74)]. 

Reweighting of the probability distribution is frequently used in MC simu- 
lations if a minus-sign problem occurs [12]. Here, the splitting into the config- 
uration weight Wh and the observable Owf is practicable provided the variance 
of both W{ and Owf is small, which is the case after the LF transformation. 
Furthermore, we require a significant overlap of the two distributions, which 
may be quantified using the Kullback-Leibler number [7], in order to avoid 
prohibitive statistical noise. In fact, our calculations show that, in general, 
for the untransformed model the reweighting method cannot be applied. For 
a detailed discussion of this point in the one-electron case see [7]. Here we 
merely note that no problems arise when simulating the transformed model. 

Apart from the significant advantage that the fermionic weight Wf only 
has to be calculated when observables are measured, the reweighting method 
becomes particularly effective in the present case when combined with the 
principal component representation introduced in section 5.4. In this case, we 
will be able to perform an exact sampling of the phonons without any auto- 
correlations. For a reliable error analysis for observables calculated according 
to equation (95) the Jackknife procedure [22] is applied. 

5.4 Principal components 

The reweighting method allows us, in principle, to skip enough sweeps be- 
tween measurements to reduce autocorrelations to a minimum. However, even 
though a single phonon update requires negligible computer time compared 
to the evaluation of Wf, for critical parameters, an enormous number of such 
steps will be necessary between successive measurements [7]. On top of that, 
reliable results require knowledge of the longest autocorrelation times, which 
have to be determined in separate simulations for each set of parameters. 

Due to the structure of the bosonic action [see equation (51)], even 
relatively small (local) changes to the phonon momenta lead to large variations 
in and hence the weight Wh- As a consequence, only minor changes may be 
proposed in order to reach a reasonable acceptance rate. Unfortunately, this 
strategy is the very origin of autocorrelations. 
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The problem can be overcome by a transformation to the normal modes of 
the phonons (along the imaginary time axis) , so that wc can sample completely 
uncorrelated configurations. As the fermion degrees of freedom are treated 
exactly, the resulting QMC method is then indeed free of any autocorrelations. 

To find such a transformation, let us recall the form of the bosonic action, 
given by equation (52), which we write as 

5b = J2pJM = EpJA'/'A^'p, (96) 

with the principal components = A^^'^p^, in terms of which the bosonic 
weight takes the simple Gaussian form 

(97) 

The sampling can now be performed directly in terms of the now variables 
^. To calculate observables we have to transform back to the physical mo- 
menta p using A~^/'^. Comparison with equation (52) shows that instead of 
the ill-conditioned matrix A we now have the ideal case that wc can easily 
generate exact samples of a Gaussian distribution. With the new coordinates 
^, the probability distribution can be sampled exactly, e.g., by the Box-Miiller 
method [23]. In contrast to a standard Markov chain MC simulation, every 
new configuration is accepted and measurements can be made at each step, 
so that simulation times are significantly reduced. 

From the definition of the principal components it is obvious that an up- 
date of a single variable ^j,,- , say, actually corresponds to a change of all pi^r' , 
t' = 1, . . . ,L. Thus, in terms of the original phonon momenta p, the updating 
becomes non-local. 

The principal component representation can be used for one, two and 
many electrons, since the bosonic action [equation (97)] is identical. This even 
holds for models including, e.g., spin-spin interactions, as long as the phonon 
operators enter in the same form as in the Holstein model. 

An important point is the combination of the principal components with 
the rewcighting method. Using the latter, the changes to the original mo- 
menta p, which are made in the simulation, do not depend in any way on the 
electronic degrees of freedom. Thus we are actually sampling a set of indepen- 
dent harmonic oscillators, as described by Sh- The crucial requirement for the 
success of this method is the use of the LF transformed model, in which the 
(bi-)polaron effects are separated from the zero-point motion of the oscillators 
around their current equilibrium positions. 

Finally, as there is no need for a warm-up phase, and owing to the sta- 
tistical independence of the configurations, the present algorithm is perfectly 
suited for parallelization. 
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Fig. 1. Average sign (sign) in the many-electron case as a function of el-ph coupling 
A in D = 1 (a) for different inverse temperatures and (b) for different values of 
the adiabaticity ratio a. Lines are guides to the eye, and errorbars are smaller than 
the symbols shown. The data presented in figures 1 and 2 are for At = 0.05. [Taken 
from [11].] 



5.5 Minus-sign problem 

The motivation for our development of a novel QMC approach to Holstein 
models was to improve on the performance of existing methods, especially 
in the many-electron case. As pointed out in [10], the LF transformation 
caiises a sign problem even for the pure Holstein model which, in general, 
may significantly affect the applicability of the method. Therefore, we briefly 
discuss the resulting limitations, focussing on the many-electron case. 

We shall sec that there is a fundamental difference between simiilations 
for one or two electrons — the carrier density being zero in the thermodynamic 
limit — and grand-canonical calculations at finite density n > 0. Whereas for 
one or two carriers the sign problem turns out to be rather uncritical — the 
average sign approaches unity upon increasing system size, in contrast to the 
usual behaviour [12] — restrictions are encountered in simulations of the many- 
electron case. 

Since Wb is strictly positive, we define the average sign as 

(sign) = («;f)b/(|w;f|)b. (98) 

For simplicity, we first show results for n = 0.5, while the effect of band 
filling will be discussed later. The choice n = 0.5 is convenient since we know 
the chemical potential, and we shall see below that the sign problem is most 
pronounced for a half-filled band. Moreover, most existing QMC results for 
the spinless Holstein model are for half filling (see references in [7]). 

Figure 1(a) shows the dependence of (sign) on the el-ph coupling strength. 
It takes on a minimum near A = 1 (for a < 1) that becomes more pronounced 
with decreasing temperature. At weak coupling (WC) and SC, (sign) !v 1, so 
that accurate simulations can be carried out. These results are quite similar 
to the cases of one or two electrons [24] . 



22 



Martin Hohenadler and Wolfgang von der Linden 




Fig. 2. Average sign in tlie many-electron case as a function of (a) band filling n, 
and (b) system size N. 

The dependence on plionon frequeney [figure 1(b)] also bears a close re- 
semblance to the polaron problem [24] . Whereas (sign) becomes very small for 
a <C 1, it increases noticeably in the non-adiabatic regime a > 1, permitting 
efficient and accurate simulations. 

As illustrated in figure 2(a), the average sign depends strongly on the band 
filling n. While it is close to one in the vicinity of n = or n = 1 (equivalent 
to one or two electrons), a significant reduction is visible near half filling 
n = 0.5. The minimum occurs at n = 0.5, and the results display particle-hole 
symmetry as expected. Here we have chosen 0t = 8, a = 0.4 and A = 1, for 
which the sign problem is most noticeable according to figure 1. 

In figure 2(b), we report the average sign as a function of system size, 
again for n = 0.5. The dependence is strikingly different from the one-electron 
case. While in the latter (sign) — !■ 1 as iV ^ oo [7,24], here the average 
sign decreases nearly exponentially with increasing system size, a behaviour 
well-known from QMC simulations of Hubbard models [12]. Obviously, this 
limits the applicability of our method. However, we shall see below that we 
can nevertheless obtain accurate results at low temperatures, small phonon 
frequencies, and over a large range of the el-ph coupling strength. Moreover, 
we would like to point out that for such parameters, other methods suffer 
strongly from autocorrelations, rendering simulations extremely difficult. 

The dependence of the sign problem on the dimension of the system is 
again similar to the single-electron case [24]. The minimum at intermediate 
A becomes more pronounced for the same parameters N, a, f3t and A as one 
increases the dimension of the cluster. 

To conclude with, wc would like to point out that, in principle, the sign 
problem can be compensated by performing sufficiently long QMC runs, but 
we have to keep in mind that the statistical errors increase proportional to 
(sign)~^ [12], setting a practical limit to the accuracy. 
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5.6 CompEirison with other approaches 

The QMC method presented above seems to be most advantageous — as com- 
pared to other approaches — in the case of the spinless Holstein model with 
many electrons. For the latter, other methods are severely restricted by auto- 
correlations, rendering accurate simulations in the physically important adi- 
abatic, IC regime virtually impossible even at moderately low temperatures. 
In contrast, the present method enables us to study the single-particle spec- 
trum on rather large clusters and for a wide range of model parameters and 
band filling (see section 6.3). Unfortunately, the generalization to the spinful 
Hubbard- Holstein model suffers severely from the sign problem. 

For the polaron and the bipolaron problem, our method requires more 
computer time than other QMC algorithms [25-28]. However, we are able to 
consider practically all parameter regimes on reasonably large clusters in one 
(polaron and bipolaron problem) and two dimensions (polaron problem). 

Finally, a discussion of the scaling of computer time with the system pa- 
rameters can be found in [7, 9, 10]. 

6 Selected results 

We now come to a selection of results obtained with the methods discussed so 
far, most of which have been published before [7,9 11]. Note that crrorbars 
will be suppressed in the figures if smaller than the symbolsize. Moreover, 
lines connecting data points are guides to the eye only. 

6.1 Small-polaron cross-over 

The Holstein model with a single electron (for a review see [29]) exhibits a 
cross-over from a large polaron (D = 1) or a qiiasi-free electron (D > 1) to a 
small polaron with increasing el-ph coupling strength. 

Qucintum Monte Carlo 

To investigate the small-polaron cross-over, following previous work [4, 16, 18- 

20,25,30 32], wc calculate the electronic kinetic energy E'kin given by equa- 
tion (75). As we shall compare results for different dimensions, we define the 
normalized quantity 

^kin = E^in/{-2tB) (99) 

with E'kin = 1 for T = and A = 0. 

The inverse temperature will be fixed to pt = 10, low enough to identify 
the cross-over. Calculations at even lower temperatures can easily be done 

for a > 1, but a < 1 requires very large numbers of measurements to en- 
sure satisfactorily small statistical errors. System sizes were 32 sites in ID, a 
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X 

Fig. 3. Normalized kinetic energy i?kin [equation (99)] of the Holstein model with 
one electron from QMC as a function of el-ph coupling A for different adiabaticity 
ratios a and different dimensions D of the lattice (A'" denotes the linear cluster size). 
Here and in subsequent figures, QMC data have been extrapolated to At = (see 
section 5.2). [Taken from [10].] 

12 X 12 cluster in 2D, and a 6 x 6 x 6 lattice in 3D. In contrast to D = 1, 2, 
where results are well converged with respect to system size, non-negligible 
finite-size effects (maximal relative changes of up to 20 % between N — 5 and 
= 6 for a <C 1; much smaller changes otherwise) are observed in three di- 
mensions. Moreover, for small N, effects due to thermal population of states 
with non-zero momentum k — absent in ground-state calculations — arc visi- 
ble, as discussed below. Nevertheless, the main characteristics are well visible 
already for N = 6. For a detailed study of finite-size and finite-temperature 
effects sec [24]. 

Figure 3 shows i?kin as a function of the el-ph coupling A for different 
phonon frequencies varying over two orders of magnitude, in one to three 

dimensions. Generally, the kinetic energy is large at WC, where the ground 
state consists of a weakly dressed electron (D > 1) or a large polaron (D = 
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1). It reduces more or less strongly — depending on a — in the SC regime, 
where a small, heavy polaron exists, defined as an electron surrounded by 
a lattice distortion essentially localized at the same site. The finite values 
of -Bkin even for large A are a result of undirected motion of the electron 
inside the surrounding phonon cloud. In contrast, the quasiparticlc weight is 
exponentially reduced in the SC regime (see, e.g., [20]), whereas the effective 
mass becomes exponentially large. 

In all dimensions, the phonon frequency has a crucial influence on the 
behaviour of the kinetic energy. While in the adiabatic regime a < 1 the 
small-polaron cross-over is determined by the condition A = E-p/2t'D > 1, 
the corresponding criterion for a > 1 is = Ep/loq > 1. The former condi- 
tion reflects the fact that the loss in kinetic energy of the electron has to be 
outweighed by a gain in potential energy in order to make small-polaron for- 
mation favourable. The latter condition expresses the increasing importance 
of the lattice energy for a > 1, since the formation of a "localized" state 
requires a sizable lattice distortion. As a consequence, for large phonon fre- 
quencies, the critical coupling shifts to Ac > 1, whereas for a < 1 we have 
Ac = 1. Additionally, the decrease of £^kin at Ac becomes signiflcantly sharper 
with decreasing phonon frequency. 

Concerning the effect of dimensionality, figure 3 reveals that, for fixed a, 
the small-polaron cross-over becomes more abrupt in higher dimensions, with 
a very sharp decrease in 3D. Nevertheless, there is no real phase transition 
[33]. Figure 3 also contains results for = 6 in one and two dimensions, i.e., 
for the same linear cluster size as in 3D (dashed lines). Clearly, for such small 
clusters, the spacing between the discrete allowed momenta k is too large 
to permit substantial thermal population, so that results are closer to the 
ground state [e.g., i?kin(A = 0) ~ 1], and exhibit a slightly more pronounced 
decrease near the critical coupling. However, the sharpening of the latter with 
increasing dimensionality is still well visible. 

Variational approach 

To test the validity of the variational approach of section 4 we have calculated 
the total energy [equation (24)] and the quasiparticlc weight [equation (29)] on 
a chistcr with = 4 for various values of a. A comparison with exact diago- 
nalization results [34] is depicted in figure 4. We only consider the regime a > 1 
where the zero-phonon approximation is expected to be justified. The overall 
agreement is strikingly good. Minor deviations from the exact results increase 
with decreasing a. For the smallest frequency shown, a= 1, the result of the 
HLF approximation is also reported. Clearly, the variational method repre- 
sents a significant improvement over the HLF approximation, underlining the 
importance of taking into account non-local distortions. Similar conclusions 
can be drawn for larger system sizes (see figure 3 in [7]). 

In figure 5 we present results for the variational displacement fields 7,5, 
which provide a measure for the polaron size. For a = 0.1 we see an abrupt 




Fig. 4. Total energy E (a) and quasiparticle weight zq (b) for A'^ = 4 as functions 
of the el-ph coupling A for different values of the adiabaticity ratio a. Symbols 

correspond to variational results and full lines represent exact T = data obtained 
with the Lanczos method [34]. Dashed lines are results of the HLF approximation. 
[Taken from [7].] 
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Fig. 5. Polaron-size parameter 74 for A'^ = 16 as a function of the el-pli coupling 
A for various distances S in the (a) adiabatic and (b) anti-adiabatic regime. Also 
shown is the LF parameter 7 [equation (10)]. [Taken from [7].] 



cross-over from a large to a small polaron at A » 1.2. For smaller A, the 
electron induces lattice distortions at neighboring sites even at a distance of 
more than three lattice constants. Above A « 1.2 we have a mobile small 
polaron extending over a single site only. In contrast, for the anti-adiabatic 
case a = 4, the cross-over is much more gradual, and 71 > even for A 1. 
The same behaviour has been found by Marsiglio [35] who determined the 
correlation function (riiXi+s) by exact diagonalization; within the variational 
approach (riiXi+s) = 15- Although in Marsiglio's results the cross-over to a 
small polaron for a = 0.1 occurs at a smaller value of the coupling A « 1, the 
simple variational approach reproduces the main characteristics. 

6.2 Bipolaron formation in the extended Holstein-Hubbard model 

In contrast to Cooper pairing of electrons with opposite momentum, two elec- 
trons may also form a bound state by travelling sufficiently close in real space. 
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Table 1. Conditions for the existence of different singlet bipolaron states in the 

onc-dimciisioiial H()lst('iii-Hul)l)ard model '9]. 
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U < 2Ep 


g<0.5 


g>0.5 


U > 4Ep (SC) 


U < 4Ep (SC) 





Bipolaron formation may be studied in the framework of the ID extended 
Holstein-Hubbard model, and a brief review of previous work has been given 
in [9,36]. Here we merely note that depending on the choice of parameters, 
the ground state of the model may either consist of two polarons, a large 
bipolaron, an inter-site bipolaron or a small bipolaron (in the singlet case). 
A summary of the conditions on the model parameters is given in table 1. 
Whereas existing work is almost exclusively concerned with the singlet case, 
here we shall also consider two electrons of the same spin. Triplet bipolarons 
are expected to play a role, e.g., in the ferromagnetic state of the manganites 
[37-39]. Furthermore, we are not aware of any previous work for y > 0. 

Quantum Monte Carlo 

Owing to the increased numerical effort compared to the onc-elcctron case, 
we shall only present results for < 12 in one dimension. However, finite-size 
effects are small even for the most critical parameters [9] . 

We define the effective kinetic energy of the two electrons as 

E^u. = E^J{-At) . (100) 

In figure 6(a) we depict E'kin as a function of the el-ph coupling for different 
values of a and U/t, at f3t = 10, i.e., much closer to the ground state than in 
some previous work [26]. 

Figure 6(a) reveals a strong decrease of i^kin near A = 0.5 for a = 0.4 
and U/t = 0. With increasing a, the cross-over becomes less pronounced, 
and shifts to larger values of A. For the same value of a, the cross-over to a 
small bipolaron is sharper than the small-polaron cross-over [cf figure 3(a)]. 
For finite on-site repulsion U/t = 4, -Ekin remains fairly large up to A w 1 (for 
a = 0.4), in agreement with the SC result Ac = 1 for [//f = 4 (see discus- 
sion in [36]). At even stronger coupling, the Hubbard repulsion is overcome, 
and a small bipolaron is formed. Again, the critical coupling increases with 
phonon frequency. Finally, the kinetic energy in the triplet case (correspond- 
ing to U/t = oo) is comparable to the results for iJ/t = 4 up to A « 1, but 
significantly larger in the SC regime since on-site bipolaron formation is not 
possible. 
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Fig. 6. Normalized kinetic energy i?kin [equation (100)] from QMC as a function 
of the el-ph coupling A for different values of the adiabaticity ratio a, the on-site 
repulsion U and the nearest-neighbour repulsion V. [(a) taken from [9].] 



The influence of nearest-neighour repulsion V is revealed in figure 6(b), 
again for U/t = 4. For all values of a shown, the cross-over sharpens noticeably 
for V > 0. The reason is that V > suppresses the (more mobile) inter-site 
bipolaron state, leading to a direct cross-over from a large to a small bipolaron. 

The nature of the bipolaron state is revealed by the correlation function 
p{S) [equation (81)], which gives the probability for the two electrons to be 
separated by a distance 5 > 0, and provides a measure of the bipolaron size. 
The phonon frequency determines the degree of retardation of the el-ph inter- 
action, and thereby limits the distance between the two electrons in a bound 
state. In the sequel, wc shall focus on the most interesting case of small phonon 
frequencies, which has often been avoided in previous work for reasons out- 
lined in section 5. 

Starting with U <C Ep , a cross-over from a small to an inter-site bipolaron 
to two weakly bound polarons takes place upon increasing the Hubbard inter- 
action [40] . Since the latter competes with the retarded el-ph interaction, the 
phonon frequency is expected to be an important parameter. In figure 7, we 
show the kinetic energy and the correlation function p{6) as a function oiU /t 
for IC A = 1. Starting from a small bipolaron for fZ/f = 0, the kinetic energy 
increases with increasing Hubbard repulsion, equivalent to a reduction of the 
effective bipolaron mass [40,41]. Although the cross-over is slightly washed 
out by the finite temperature in our simulations, there is a well- conceivable 
increase in i?kin up to U /t 4, above which the kinetic energy begins to 
decrease slowly. The increase of i?kin originates from the breakup of the small 
bipolaron, as indicated by the decrease of p(0) in figure 7(b). Close to U/t = 4, 
the curves for p(0) and p{l) cross, and it becomes more favourable for the two 
electrons to reside on neighboring sites. The inter-site bipolaron only exists 
below a critical Hubbard repulsion Uc- The latter is given by Uc = 2Ep (i.e., 
here Uc/t = 4) at weak el-ph coupling, and by Uc = 4i?p at SC. For an inter- 
mediate value A = 1 as in figure 7, the cross-over from the inter-site state to 
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Fig. 7. (a) Normalized kinetic energy iSkin and (b) correlation functions p(0), p(l) 
from QMC as a function of the Hubbard repulsion U jt for different values of the 
adiabaticity ratio a. [Taken from [9].] 



two weakly bound polarons is expected to occur somewhere in between, but 
is difficult to locate exactly from the QMC results. 

Figure 7 further illustrates that the cross-over becomes steeper with de- 
creasing phonon frequency. In the adiabatic limit a = 0, it has been shown 
to be a first-order phase transition [42], whereas for a > retardation effects 
suppress any non-analytic behaviour. At the same J7/t, i?kin increases with a 
since for a fixed A, the bipolaron becomes more weakly bound. For the same 
reason, the cross-over to an inter-site bipolaron — showing up in figure 7 as a 
crossing of /9(0) and p(l) - shifts to smaller values of U jt. 

Let us now consider the effect of temperature on p(S). To this end, we plot 
in figures 8(a) - (c) p{S) at different temperatures, for parameters correspond- 
ing to the three regimes of a large, small and inter-site bipolaron, respectively. 

For the parameters in figure 8(a) ([//f = 0, A = 0.25), the two electrons 
are most likely to occupy the same site, but the bipolaron extends over a 
distance of several lattice constants. Clearly, in this regime, the cluster size 

= 12 used here is not completely satisfactory, but still provides a fairly 
accurate description as can be deduced from calculations for A' = 14 (not 
shown). Nevertheless, on such a small cluster, no clear distinction between 
an extended bipolaron and two weakly bound polarons can be made. As the 
temperature increases from (it = 10 to /?t = 1, the probability distribution 
broadens noticeably, i.e., it becomes more likely for the two electrons to be fur- 
ther apart. In particular, for the highest temperature shown, /9(0) has reduced 
by about 30 % compared to /3f = 10. 

A different behaviour is observed for the small bipolaron, which exist at 
stronger el-ph coupling A = 1. Figure 8(b) reveals that p{6) peaks strongly 
at ^ = 0, but is very small for ^ > at low temperatures. Increasing tem- 
perature, p(6) remains virtually unchanged up to /?t = 3. Only at very high 
temperatures there occurs a noticeable transfer of probability from 5 = to 
5 > 0. At the highest temperature shown, = 0.5, the two electrons have a 
non-negligible probability for traveling a finite distance 5 > apart, although 
most of the probability is still contained in the peak located at (5 = 0. 
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Fig. 8. Correlation function p{S) from QMC as a function of S for different inverse 
temperatures (5, N = 12 and a = 0.4. [Taken from [9].] 

Finally, we consider in figure 8(c) the inter-site bipolaron, taking U /t = 4 
and A = 1 (of figure 2 in [9]). At low temperatures, p{6) takes on a maximum 
for 5 = 1. For smaller values of the latter diminishes, until at (3t = 1, the 
distribution is completely fiat, so that all 5 are equally likely. 

The different sensitivity of the bipolaron states to changes in temperature 
found above can be explained by their different binding energies. The latter 
is given by AEq, = E^P — 2E[^\_ where e'"^^ and e'"^^ denote the ground-state 
energy of the model with one and two electrons, respectively. 

Generally, the thermal dissociation is expected to occur at a temperature 
such that the thermal energy k^T = (jST)^^ becomes comparable to AEq, in 
accordance with our numerical data. The large and the inter-site bipolaron 
are relatively weakly bound as a result of the rather small effective interac- 
tion Ucs ~ U - 2Ep [36]. The binding energies are AEq « -(0.32 ± 0.08)< 
and —(0.28 ± 0.08)t, respectively, so that we expect a critical inverse tem- 
perature f3t w 2.5-5, in agreement with figures 8(a) and (c). In contrast. 
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Fig. 9. Variational results for the normalized kinetic energy iSkin as a function of 
the el-ph coupling A, and for different adiabaticity ratios a. Also shown are results 
of the HLF approximation. [Taken from [9].] 

the small bipolaron in figure 8(b) has a significantly larger binding energy 
AE ?a -(3.43 ± 0.09)i, and therefore remains stable up to f3t « 0.3. Ther- 
mal dissociation of bipolarons occurs at even lower temperatures for V > 0, 
especially in the triplet case, owing to the reduced binding energy. 

Variational approach 

Whereas the QMC approach is limited to finite temperatures and relatively 
small clusters, the variational method of section 4 yields ground-state results 
on much larger systems. To scrutinize the quality of the variational method, we 
compare the ground-state energy for U/t = to the most accurate approach 
currently available in one dimension, namely the variational diagonalization 
[40]. We find a good agreement over the whole range of A. As expected from 
the nature of the approximation, slight deviations occur for a < 1, similar to 
the one-electron case. 

Despite the success in calculating the total energy — being the quantity 
that is optimized — one has to be careful not to overestimate the validity of any 
variational method. To reveal the shortcomings of the current approach, we 
show in figure 9 the normalized kinetic energy Ekin = ieff [see equations (31) 
and (100)] as a function of el-ph coupling, and for different a. We have chosen 
A'' = 25 to ensure negligible finite-size effects. In principle, figure 9 displays 
a behaviour similar to the QMC data in figure 6(a). There is a jump-like 
decrease of E'kin near A = 0.5 for a = 0.4, which becomes washed out and 
moves to larger A with increasing plionon frequency. For a ~ 0.4, the cross- 
over in the variational results is much too steep, regardless of the fact that 
the latter are for T = 0, a common defect of variational methods. Moreover, 
for a = 0.4 2, the variational kinetic energy is too small above the bipolaron 
cross-over compared to the QMC data, whereas for a = 4, the decay of i^kin 
with increasing A is too slow. 

The reason for the failure is the absence of retardation effects, which play 
a dominant role in the formation of bipolaron states. The increased impor- 
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Fig. 10. One-electron spectral function A{k, u> — from QMC for different band 
fillings n, N = 32, /3t = 8, a = 0.4, and A = 0.1. Here and in subsequent figures 
Zir = 0.1. [Taken from [11].] 



tance of the phonon dynamics — not included in the variational method — for 
the two-electron problem leads to a less good agreement with exact results 
than in the one-electron case. In particular, our variational results overesti- 
mate the position of the cross-over (figure 9) compared to the value Ac = 0.5 
expected in the adiabatic regime. Nevertheless, the method represents a signif- 
icant improvement over the simple HLF approximation, due to the variational 
determination of the parameters 7^ . This is illustrated in figure 9, where we 
also show the HLF result E^in = for a = 0.4 and 4.0. In contrast to 
the variational approach, the HLF approximation yields an exponentially de- 
creasing kinetic energy for all values of the phonon frequency. Whereas such 
behaviour actually occurs in the anti- adiabatic limit a — > 00, the situation 
is different for small a [sec figures 6(a) and 9]. The variational method pre- 
sented here accounts qualitatively for the influence of the phonon frequency 
on bipolaron formation. 



6.3 Many-polciron problem 

We review recent results on the carrier-density dependence of photocmission 
spectra of many-polaron systems in the framework of the spinless Holstein 
model (63) in one dimension. We shall see that the sensitivity to changes 
in n strongly depends on the phonon frequency and el-ph coupling strength, 
with the most interesting physics being observed in the adiabatic, IC regime 
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Fig. 11. One-electron density of states N{uj — jj.) from QMC for different band 
fillings n, N = 32, pt = 8, a = 0.4 and A = 0.1. [Taken from [11].] 



often realized experimentally. This regime is characterized by the existence of 
large polarons at low carrier density. At larger densities, a substantial overlap 
of the single-particle wavefunctions occurs, leading to a dissociation of the 
individual polarons and finally to a restructuring of the whole many-particle 
ground state. Note that the many-polaron problem has since been studied also 
by means of other methods [43-45], confirming the original findings of [11]. 

Weak coupling 

For WC A = 0.1, the sign problem is not severe (section 5.5) so that simu- 
lations can easily be performed for large lattices with A'^ = 32, making the 
dispersion of quasiparticle features well visible. 

Figure 10 shows the evolution of the one-electron spectral function A{k, u— 
fi) with increasing electron density n. At first sight, we see that the spectra 
bear a close resemblance to the free-electron case, i.e., there is a strongly 
dispersive band running from —2t to 2t which can be attributed to weakly 
dressed electrons. As expected, the height (width) of the peaks increases (de- 
creases) significantly in the vicinity of the Fermi momentum fcp, determined 
by the crossing of the band with the chemical potential. However, in contrast 
to the case of a rigid tight- binding band, we shall see below (figure 11) that a 
significant redistribution of spectral weight occurs with increasing n. 

We would like to point out that the apparent absence of any phonon signa- 
tures in figure 10 is not a defect of the maximum entropy method, but results 
from the large scale of the 2:-axis chosen. As a consequence, the peaks running 
close to the bare band dominate the spectra and suppress any small phonon 
peaks present. At higher resolution, for all densities n = 0.1-0.4, we observe 
the band flattening [46-48] at large wavevectors which originates from the in- 
tersection of the approximately free-electron dispersion with the bare phonon 
energy at w — fj, = u)q. 
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Fig. 12. One-electron spectral function A{k, u> — from QMC for different band 
fillings n, AT = 32, /3t = 8, a = 0.4, and A = 2. [Taken from [11].] 



To complete our discussion of the WC regime, we show in figure 1 1 the one- 
electron density of states (DOS) N{uj — ji) given by equation (94). Clearly, 
for small n, there is a peak with large spectral weight at the Fermi level. 
In contrast, for large n, the tendency toward formation of a Peierls- (band-) 
insulating state at n = 0.5 suppresses the DOS at the Fermi level, although we 
are well below the critic;al value of A at which the cross-over to the insulating 
state takes place at T = [49,50]. The additional small features separated 
from 12 by the bare phonon energy ivq will be discussed below. 

Strong coupling 

We now turn to the SC limit taking A = 2. At low density n = 0.1 [fig- 
ure 12(a)], we expect the well-known, almost flat polaron band having ex- 
ponentially reduced spectral weight (given by in the single-electron, SC 
limit) which, nevertheless, can give rise to coherent transport at T = 0. As 
discussed in [11], such weak signatures are difficult to determine accurately 
using the maximum entropy method. Generally, it is known that the reliability 
of dynamic properties obtained by means of the maximum entropy method 
crucially depends on the size of statistical errors and the general structure of 
the spectra. A detailed discussion of this point has been given in [11]. 

Besides, the spectrum consists of two incoherent features located above and 
below the chemical potential, which reflect the phonon-mediated transitions to 
high-energy electron states. Here, the maximum of the photoemission spectra 
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Fig. 13. One-electron density of states A'^(a; — /i) from QMC for different band 
filhngs n and cluster sizes N, f3t = 8, a = 0.4 and A = 2. [Taken from [11].] 



(w — /X > 0) follows a tight-binding cosine dispersion. The incoherent part of 
the spectra is broadened according to the phonon distribution. 

For all band fillings, the chemical potential is expected to be located in 
a narrow polaron band with little spectral weight. There exists a finite gap 
to the photoemission (inverse photoemission) parts of the spectrum, so that 
the system typifies as a polaronic metal. We shall see below that a completely 
different behaviour is obscrvc;d at IC. Notice that the incoherent inverse pho- 
toemission (photoemission) signatures are more pronounced at small (large) 
wavevectors. 

Finally, for n = 0.4 [figure 12(d)], the incoherent features lie rather close 
to the Fermi level, thus being accessible by low-energy excitations. Now, the 
photoemission spectrum for < 7r/2 is almost symmetric to the inverse pho- 
toemission spectrum for k > n/2 and already reveals the gapped structure 
which occurs at n = 0.5 due to charge- density- wave formation accompanied 
by a Peierls distortion [50]. 

As in the WC case discussed above, the properties of the system also 
manifest itself in the DOS, shown in figure 13. Owing to the strong el-ph 
interaction, the spectral weight at the chemical potential is exponentially small 
for all fillings n. At half filling, the DOS exhibits particle-hole synunctry, and 
the system can be described as a Peierls insulator, consisting of a polaronic 
superlattice. In contrast to the WC case, the ground state is characterized as 
a polaronic insulator rather than as a band insulator. 

Intermediate coupling 

As discussed in the introduction, a cross-over from a polaronic state to a sys- 
tem with weakly dressed electrons can be expected in the IC regime. Here we 
choose A = 1, which corresponds to the critical value for the small-polaron 
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Fig. 14. One-electron spectral function A{k, u> — from QMC for different band 
fillings n and cluster sizes N, f3t = 8, a = 0.4, and A = 1. [Taken from [11].] 



cross-over in the one-electron problem [cf figure 3(a)]. Owing to the sign prob- 
lem, which is particularly noticeable for A = 1 (see figure 1), we have to 
decrease the system size as we increase the electron density n. 

We shall see that the cross-over is rather difficult to detect from the QMC 
results only. However, the data presented here are perfectly consistent with 
more recent studies employing other methods such as exact diagonalization 
[11], cluster perturbation theory [44] or self-energy calculations [43]. 

Figure 14 shows the spectral function for A = 1 and increasing band filling. 
Owing to the overlap of large polarons in the IC regime, we start with a very 
low density n = 0.05 [figure 14(a)]. Compared to the behaviour for A = 
2 [figure 12(a)], we notice that the polaron band now lies much closer to 
the incoherent features, and that there is a mixing of these two parts of the 
spectrum at small values of k. Nevertheless, the almost flat polaron band is 
well visible for large k. 

With increasing density, the polaron band merges with the incoherent 
peaks at higher energies, signaling the above-anticipated density-driven cross- 
over from a polaronic to a (diffusive) metallic state, with the broad main band 
crossing the Fermi level. 

Further information about the density dependence can be obtained from 
the one-electron DOS. The latter is presented in figure 15 for different fillings 
n = 0.05-0.5. As in figure 14, the cluster size is reduced with increasing n in 
order to cope with the sign problem. To illustrate the rather small influence 
of finite-size effects, figure 15 also contains results for A'' = 10. 
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Fig. 15. One-electron density of states N{lu — ji) from QMC for different band 
fillings n, cluster sizes N and inverse temperatures (5. Here a = 0.4 and A = 1. 
[Taken from [11].] 

For low density n = 0.05, the DOS in figure 15 lies in between the results 
for WC and SC discussed above. Although the spectral weight at the chemical 
potential is strongly reduced compared to A = 0.1, N{0) is still significantly 
larger than for A = 2. 

When the density is increased to n = 0.2, the DOS at the chemical poten- 
tial increases, as a result of the dissociation of polarons. Increasing n further, a 
pseudogap begins to form at /.t, which is a precursor of the charge-density-wave 
gap at half filling and zero temperature. 

In the case of half filling n — 0.5, the DOS has become symmetric with 
respect to fi. There are broad features located either side of the chemical 
potential, which take on maxima close to u — /i = ±Ep. However, apart from 
the SC case, where the single-polaron binding energy is still a relevant energy 
scale, the position of these peaks is rather determined by the energy of the 
upper and lower bands, split by the formation of a Peierls state. The gap of 
size ^ A expected for the insulating charge-ordered state at T = is partially 
filled in due to the finite temperature considered here. 

Furthermore, we find additional, much smaller features roughly separated 
from /i by the bare phonon frequency loq, whose height decreases with de- 
creasing temperature, as revealed by the results for (3t = 10 (figure 15). These 
peaks — not present at T = [50,51] — arise from thermally activated tran- 
sitions to states with additional phonons excited, and are also visible in fig- 
ures 11 and 13. While for WC (A = 0.1, figure 11), the maximum of these 
features is almost exactly located at [a; — — loq, it moves to [w — « 1.25aJo 
for IC (A = 1, figure 15), and finally to \uj - n\ ^ 2.5wo for SC (A = 2, fig- 
ure 13). Although the exact positions of the peaks are subject to uncertainties 
due to the maximum entropy method, this evolution refiects the shift of the 
maximum in the phonon distribution function with increasing coupling. The 
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maximum entropy method yields an envelope of the multiple peaks separated 
by Wq. 

Anti-adiabatic regime 

The comparison of the spectral functions for n = 0.1 and n = 0.3 in figure 10 of 
[11] reveals that there is no density- driven cross-over of the system as observed 
in the adiabatic case even for the critical value = 1. In particular, owing 
to the large phonon energy, there are no low-energy excitations close to the 
polaron band, so that the latter remains well separated from the incoherent 
features even for n = 0.3. Furthermore, the spectral weight of the polaron 
band also remains almost unchanged as we increase the density from n = 0.1 
to n = 0.3. Consequently, almost independent small polarons are formed also 
at finite electron densities, in accordance with previous findings for small 
systems [52]. 

7 Summary 

We have reviewed quantum Monte Carlo and variational approaches to Hol- 
stein models based on Lang-Firsov transformations of the Hamiltonian. The 
methods have been applied to investigate single polarons and bipolarons, re- 
spectively, as well as a many-polaron system. 

The variational methods include displacements of the lattice at all lat- 
tice sites, which enables them to quite accurately describe large polaron or 
bipolaron states. 

Using the transformed Hamiltonian, we have shown that quantum Monte 
Carlo simulation can be based on exact sampling without autocorrelations, 
which proves to be an enormous advantage for small phonon frequencies or 
low temperatures. Indeed, we have used a grand-canonical algorithm to obtain 
dynamical properties of many-polaron systems in all interesting parameter 
regimes. Such simulations are currently not possible with other Monte Carlo 
methods. 
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